"""
CBMPy: Constraint Based Modelling in Python (http://pysces.sourceforge.net/cbm)
============
Copyright (C) 2009-2012 Brett G. Olivier, VU University Amsterdam, Amsterdam, The Netherlands

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>

Author: Brett G. Olivier
Contact email: bgoli@users.sourceforge.net
Last edit: $Author: bgoli $ ($Id: CBWrite.py 990 2012-02-15 13:52:11Z bgoli $)

"""

if vars().has_key('cDir'):
    cDir=vars()['cDir']
else:
    cDir=None
import os, time, cStringIO, numpy
import CBTools, CBXML
from CBVersion import __DEBUG__, __version__

_HAVE_SYMPY_ = None
try:
    import sympy
    _HAVE_SYMPY_ = True
except ImportError:
    print 'Rational IO not available'
    _HAVE_SYMPY_ = False

##  _PREC_ = '%f'

def writeSBML3(fba, fname, directory=None, sbml_level_version=None):
    """
    Takes an FBA model object and writes it to file as SBML L2 with FBA annotations:

     - *fba* an fba model object
     - *fname* the model will be written as XML to *fname*
     - *sbml_level_version* [default=None] a tuple containing the SBML level and version e.g. (2,4)

    This is a utility wrapper for the function `CBXML.sbml_writeSBML2FBA`

    """
    print '\n*****\nThis function is being deprecated, please use writeSBML2FBA() instead!\n*****\n'
    CBXML.sbml_writeSBML2FBA(fba, fname, directory, sbml_level_version)

def writeSBML2FBA(fba, fname, directory=None, sbml_level_version=None):
    """
    Takes an FBA model object and writes it to file as SBML L2 with FBA annotations.
    Note if you want to write BiGG/FAME style annotations then you must use *sbml_level_version=(2,1)*

     - *fba* an fba model object
     - *fname* the model will be written as XML to *fname*
     - *sbml_level_version* [default=None] a tuple containing the SBML level and version e.g. (2,4)

    This is a utility wrapper for the function `CBXML.sbml_writeSBML2FBA`

    """
    CBXML.sbml_writeSBML2FBA(fba, fname, directory, sbml_level_version)


def writeSensitivitiesToCSV(sensitivities, fname):
    """
    Write out a sensitivity report using the objective sensitivities and
    bound sensitivity dictionaries created by e.g. cplx_getSensitivities().

     - *sensitivity* tuple containing

      - *obj_sens* dictionary of objective coefficient sensitivities (per flux)
      - *rhs_sens* dictionary of constraint rhs sensitivities (per constraint)
      - *bound_sens* dictionary of bound sensitivities (per flux)

     - *fname* output filename e.g. fname.csv

    """
    obj_sens = sensitivities[0]
    rhs_sens = sensitivities[1]
    bound_sens = sensitivities[2]
    F = file(fname+'_flux_sensitivity.csv', 'w')
    head = "Flux,Reduced cost,OCS low,OC value,OCS high,LB low,LB high,UB low,UB high"
    F.write(head+'\n')
    for j in obj_sens.keys():
        rc = obj_sens[j][0]
        lcs = obj_sens[j][1]
        ocv = obj_sens[j][2]
        ucs = obj_sens[j][3]
        if j in bound_sens:
            lbs = bound_sens[j][0]
            lb = bound_sens[j][1]
            ub = bound_sens[j][2]
            ubs = bound_sens[j][3]
        else:
            lbs = 0
            lb = 0
            ub = 0
            ubs = 0
        F.write('%s,%s,%s,%s,%s,%s,%s,%s,%s\n' % (j,rc,lcs,ocv,ucs,lbs,lb,ub,ubs))

    for j in bound_sens.keys():
        if j not in obj_sens:
            rc = 0
            lcs = 0
            ocv = 0
            ucs = 0
            lbs = bound_sens[j][0]
            lb = bound_sens[j][1]
            ub = bound_sens[j][2]
            ubs = bound_sens[j][3]
            F.write('%s,%s,%s,%s,%s,%s,%s,%s,%s\n' % (j,rc,lcs,ocv,ucs,lbs,lb,ub,ubs))
    F.flush()
    F.close()

    F = file(fname+'_constraint_sensitivity.csv', 'w')
    F.write('Constraint,RHS low,RHS,RHS high\n')
    for c in rhs_sens.keys():
        F.write('%s,%s,%s,%s\n' % (c, rhs_sens[c][0], rhs_sens[c][1], rhs_sens[c][2]))
    F.flush()
    F.close()

def WriteModelRaw(fba, work_dir=cDir):
    """
    Writes a fba (actually just dumps it) to a text file.

     - *fba* an instantiated FBAmodel instance
     - *work_dir* directory designated for output

    """
    if work_dir == None:
        work_dir = os.getcwd()
    FF = file(os.path.join(work_dir,'WriteModelRawOutput.txt'), 'w')
    FF.write('Species information:\n\n')
    for s in fba.species:
        FF.write('%s: value=%f, is_boundary=%s, name=%s\n' % (s.getPid(), s.value, int(s.is_boundary), s.name))
    FF.write('\nReaction information:\n\n')
    for r in fba.reactions:
        FF.write('%s: reversible=%s, name=%s\n' % (r.getPid(), int(r.reversible), r.name))
        FF.write('\t%s\n' % r.getReagentRefs())
    FF.write('\nConstraint information:\n\n')
    for c in fba.flux_bounds:
        FF.write('%s %s %f\n' % (c.reaction, c.operation, c.value))
    FF.write('\nObjective information:\n\n')
    FF.write('Active Objective: %s (%s)\n' % (fba.objectives[fba.activeObjIdx].getPid(), fba.objectives[fba.activeObjIdx].operation))
    for o in fba.objectives:
        FF.write('%s: %s, %s\n' % (o.getPid(), o.operation, o.getFluxObjectiveReactions()))
    if hasattr(fba.N, 'shape'):
        FF.write('\nStoichiometric information:\n\n')
        FF.write('N-matrix dimensions = (%s,%s)\n' % fba.N.shape)
    FF.close()
    print 'WriteModelRaw has written a file to %s' % os.path.join(work_dir,'WriteModelRawOutput.txt')

def BuildLPFluxBounds(fba, use_rational=False):
    """
    Build and return a cStringIO that contains the flux bounds in LP format
    """
    if not _HAVE_SYMPY_ and use_rational:
        use_rational = False
        print 'Warning: install Sympy for rational IO'

    FFS = cStringIO.StringIO()
    c2s = {}
    for c in fba.flux_bounds:
        if __DEBUG__: print '%s: %s' % (c.getPid(), c.operation)
        minv = None
        maxv = None
        oper = None
        R = c.reaction
        if c.operation in ['less','lessEqual']:
            maxv = c.value
            oper = c.operation
            if oper == 'less':
                oper = '<='
                if __DEBUG__: print 'LP Bounds operator must be <= not <'
            else:
                oper = '<='
        if c.operation in ['greater','greaterEqual']:
            minv = c.value
            oper = c.operation
            if oper == 'greater':
                oper = '<='
                if __DEBUG__: print 'LP Bounds operator must be <= not <'
            else:
                oper = '<='
        if c.operation == 'equals':
            maxv = c.value
            oper = '='
        if maxv != None:
            if use_rational:
                c2s.update({R+'r' : '%s %s' % (oper, sympy.Rational(maxv))})
            else:
                c2s.update({R+'r' : '%s %s' % (oper, maxv)})
        elif minv != None:
            if use_rational:
                c2s.update({R+'l' : '%s %s' % (sympy.Rational(minv), oper)})
            else:
                c2s.update({R+'l' : '%s %s' % (minv, oper)})

    c2sk = c2s.keys()
    c2sk.sort()
    for r in (c2sk):
        rht = None
        lht = None
        R = r[:-1]
        if c2s.has_key(R+'r'):
            rht = c2s[R+'r']
            c2s.pop(R+'r')
        if c2s.has_key(R+'l'):
            lht = c2s[R+'l']
            c2s.pop(R+'l')
        if lht != None and rht != None:
            if __DEBUG__: print '%s %s %s' % (lht, R, rht)
            ##  FFS.write('%s: %s %s %s\n' % (R, lht, R, rht))
            FFS.write('%s %s %s\n' % (lht, R, rht)) # THIS MUST BE SO for GLPK
        elif lht == None and rht == None:
            if __DEBUG__: print 'Skipping: (%s, %s, %s)' % (lht, R, rht)
        elif lht != None:
            if __DEBUG__: print '%s %s' % (lht, R)
            #FFS.write('%s: %s %s\n' % (R, lht, R))
            FFS.write('%s %s\n' % (lht, R))
        elif rht != None:
            if __DEBUG__: print '%s %s' % (R, rht)
            #FFS.write('%s: %s %s\n' % (R, R, rht))
            FFS.write('%s %s\n' % (R, rht))
        else:
            print 'CONFUSION: (%s, %s, %s)' % (lht, R, rht)

    return FFS

def BuildLPConstraints(fba, use_rational=False):
    """
    Build and return a cStringIO that contains constraint constructed from
    the StoichiometeryLP object
    
     - *fba* an fba model object which has a stoichiometry
     - *use_rational* write rational number output [default=False]
    """
    
    if not _HAVE_SYMPY_ and use_rational:
        use_rational = False
        print 'Warning: install Sympy for rational output'
    rebuild_stoich = False
    if not hasattr(fba,'N') or fba.N == None:
        rebuild_stoich = True
    else:
        if len([s for s in fba.species if not s.is_boundary]) != fba.N.array.shape[0]:
            rebuild_stoich = True
        elif len(fba.reactions) != fba.N.array.shape[1]:
            rebuild_stoich = True
    if rebuild_stoich:
        print '\nWarning FBA object has inconsistant stoichiometric matrix, rebuilding it now.'
        CBTools.addStoichToFBAModel(fba)

    constr = {}
    for r in range(fba.N.array.shape[0]):
        rowName = fba.N.row[r]
        RCon = []
        for c in range(fba.N.array.shape[1]):
            colName = fba.N.col[c]
            colCoef = fba.N.array[r,c]
            if colCoef != 0.0:
                RCon.append((colCoef, colName))
        constr.update({rowName : RCon})

    FFS = cStringIO.StringIO()
    constrsk = constr.keys()
    constrsk.sort()
    ##  for r in constrsk:
    for r in range(len(fba.N.row)):
        if len(constr[fba.N.row[r]]) > 0:
            FFS.write(' %s: ' % fba.N.row[r])
            ##  FFS.write(' ')
            for col in constr[fba.N.row[r]]:
                if use_rational:
                    if col[0] > 0.0:
                        FFS.write('+%s %s ' % (sympy.Rational('%s' % col[0]), col[1]))
                    else:
                        FFS.write('%s %s ' % (sympy.Rational('%s' % col[0]), col[1]))
                else:
                    if col[0] > 0.0:
                        FFS.write('+%s %s ' % (col[0], col[1]))
                    else:
                        FFS.write('%s %s ' % (col[0], col[1]))
            FFS.write('%s ' % fba.N.operators[r])
            FFS.write('%s\n' % fba.N.RHS[r])
    return FFS

def BuildLPUserConstraints(fba, use_rational=False):
    """
    Build and return a cStringIO that contains constraint constructed from
    the StoichiometeryLP object
    
     - *fba* an fba model object which has a stoichiometry
     - *use_rational* write rational number output [default=False]
    """
    
    if not _HAVE_SYMPY_ and use_rational:
        use_rational = False
        print 'Warning: install Sympy for rational output'
    rebuild_stoich = False
    assert fba.user_constraints != None, "\nNo user constraints to build"
    if not hasattr(fba,'CM') or fba.CM == None:
        rebuild_stoich = True

    if rebuild_stoich:
        print '\nWarning FBA object has inconsistant user constraint matrix, rebuilding it now.'
        CBTools.addStoichToFBAModel(fba)

    constr = {}
    for r in range(fba.CM.array.shape[0]):
        rowName = fba.CM.row[r]
        RCon = []
        for c in range(fba.CM.array.shape[1]):
            colName = fba.CM.col[c]
            colCoef = fba.CM.array[r,c]
            if colCoef != 0.0:
                RCon.append((colCoef, colName))
        constr.update({rowName : RCon})

    FFS = cStringIO.StringIO()
    constrsk = constr.keys()
    constrsk.sort()
    ##  for r in constrsk:
    for r in range(len(fba.CM.row)):
        if len(constr[fba.CM.row[r]]) > 0:
            FFS.write(' %s: ' % fba.CM.row[r])
            ##  FFS.write(' ')
            for col in constr[fba.CM.row[r]]:
                if use_rational:
                    if col[0] > 0.0:
                        FFS.write('+%s %s ' % (sympy.Rational('%s' % col[0]), col[1]))
                    else:
                        FFS.write('%s %s ' % (sympy.Rational('%s' % col[0]), col[1]))
                else:
                    if col[0] > 0.0:
                        FFS.write('+%s %s ' % (col[0], col[1]))
                    else:
                        FFS.write('%s %s ' % (col[0], col[1]))
            FFS.write('%s ' % fba.CM.operators[r])
            FFS.write('%s\n' % fba.CM.RHS[r])
    return FFS


def BuildLPConstraintsRelaxed(fba):
    """
    Build and return a cStringIO that contains the constaints in LP format
    Relaxed refers to dS/dt >= 0
    """
    raise DeprecationWarning, "\nThis method is deprecated"
    if not hasattr(fba,'N') or fba.N == None:
        print '\nWarning FBA object has no stoichiometric matrix constructing it now.'
        CBTools.addStoichToFBAModel(fba)
        time.sleep(1)    
    
    constr = {}
    for r in range(fba.N.array.shape[0]):
        rowName = fba.N.row[r]
        RCon = []
        for c in range(fba.N.array.shape[1]):
            colName = fba.N.col[c]
            colCoef = fba.N.array[r,c]
            if colCoef != 0.0:
                RCon.append((colCoef, colName))
        constr.update({rowName : RCon})

    FFS = cStringIO.StringIO()
    constrsk = constr.keys()
    constrsk.sort()
    for r in fba.N.row:
        if len(constr[r]) > 0:
            FFS.write('%s: ' % r)
            ##  FFS.write(' ')
            for col in constr[r]:
                if col[0] > 0.0:
                    FFS.write('+%s %s ' % (col[0], col[1]))
                else:
                    FFS.write('%s %s ' % (col[0], col[1]))
            FFS.write('>= 0\n')

    return FFS

def BuildLPConstraintsStrict(fba, use_rational=False):
    """
    Build and return a cStringIO that contains the constaints in LP format
    Strict refers to dS/dt = 0
    """
    ##  print "Consider using the new BuildLPConstraints() method"
    
    if not _HAVE_SYMPY_ and use_rational:
        use_rational = False
        print 'Warning: install Sympy for rational IO'

    if not hasattr(fba,'N') or fba.N == None:
        print '\nWarning FBA object has no stoichiometric matrix constructing it now.'
        CBTools.addStoichToFBAModel(fba)
        time.sleep(1)

    constr = {}
    for r in range(fba.N.array.shape[0]):
        rowName = fba.N.row[r]
        RCon = []
        for c in range(fba.N.array.shape[1]):
            colName = fba.N.col[c]
            colCoef = fba.N.array[r,c]
            if colCoef != 0.0:
                RCon.append((colCoef, colName))
        constr.update({rowName : RCon})

    FFS = cStringIO.StringIO()
    constrsk = constr.keys()
    constrsk.sort()
    ##  for r in constrsk:
    for r in fba.N.row:
        if len(constr[r]) > 0:
            FFS.write(' %s: ' % r)
            ##  FFS.write(' ')
            for col in constr[r]:
                if use_rational:
                    if col[0] > 0.0:
                        FFS.write('+%s %s ' % (sympy.Rational('%s' % col[0]), col[1]))
                    else:
                        FFS.write('%s %s ' % (sympy.Rational('%s' % col[0]), col[1]))
                else:
                    if col[0] > 0.0:
                        FFS.write('+%s %s ' % (col[0], col[1]))
                    else:
                        FFS.write('%s %s ' % (col[0], col[1]))
            if len(constr[r]) != 0:
                FFS.write('= 0\n')
            else:
                FFS.write('\n')

    return FFS

def BuildLPConstraintsMath(fba, use_rational=False):
    """
    Build and return a cStringIO that contains the constaints in LP format
    Strict refers to dS/dt => 0 and dS/dt <= 0
    """
    raise DeprecationWarning, "\nThis method is deprecated"
    if not _HAVE_SYMPY_ and use_rational:
        use_rational = False
        print 'Warning: install Sympy for rational IO'

    if not hasattr(fba,'N') or fba.N == None:
        print '\nWarning FBA object has no stoichiometric matrix constructing it now.'
        CBTools.addStoichToFBAModel(fba)
        time.sleep(1)

    constr = {}
    for r in range(fba.N.array.shape[0]):
        rowName = fba.N.row[r]
        RCon = []
        for c in range(fba.N.array.shape[1]):
            colName = fba.N.col[c]
            colCoef = fba.N.array[r,c]
            if colCoef != 0.0:
                RCon.append((colCoef, colName))
        constr.update({rowName : RCon})

    FFS = cStringIO.StringIO()
    constrsk = constr.keys()
    constrsk.sort()
    for r in fba.N.row:
        if len(constr[r]) > 0:
            FFS.write('%sn1: ' % r)
            ##  FFS.write(' ')
            for col in constr[r]:
                if col[0] > 0.0:
                    FFS.write('+%s %s ' % (col[0], col[1]))
                else:
                    FFS.write('%s %s ' % (col[0], col[1]))
            # This is a fudge
            FFS.write('>= 0\n')

    for r in fba.N.row:
        if len(constr[r]) > 0:
            FFS.write('%sn2: ' % r)
            for col in constr[r]:
                if -col[0] > 0.0:
                    FFS.write('+%s %s ' % (-col[0], col[1]))
                else:
                    FFS.write('%s %s ' % (-col[0], col[1]))
            # This is a fudge
            FFS.write('>= 0\n')
    return FFS

def WriteModelLPOld(fba, work_dir=None, multisymb=' ', lpt=True, constraint_mode='strict', use_rational=False, format='%s'):
    """
    Writes a fba as an LP/LPT

     - *fba* an instantiated FBAmodel instance
     - *work_dir* directory designated for output
     - *multisymb* the multiplication symbol (default: <space>)
     - *lpt* the file format (default: True for lpt) or False for lp

    """
    
    print "\nTHIS FUNCTION IS DEPRECATED\n"
    
    if not _HAVE_SYMPY_ and use_rational:
        use_rational = False
        print 'Warning: install Sympy for rational IO'

    FNAME = None
    if work_dir == None:
        FnameTmp = fba.getPid()
    else:
        FnameTmp = os.path.join(work_dir, fba.getPid())
    if use_rational:
        FnameTmp = FnameTmp+'.rat'
    if not lpt:
        FNAME = FnameTmp+'.lp'
        FF = file(FNAME, 'w')
        FF.write('Problem\n %s\n\n' % FnameTmp)
    else:
        FNAME = FnameTmp+'.lp'
        FF = file(FNAME, 'w')
        FF.write('\\\\ %s \n\n' % FnameTmp)
    objO = fba.objectives[fba.activeObjIdx].operation.lower()
    objO = objO[0].upper() + objO[1:]
    FF.write('%s\n' % objO)
    objStr = '%s_objf: ' % fba.objectives[fba.activeObjIdx].getPid()
    for fObj in fba.objectives[fba.activeObjIdx].fluxObjectives:
        sign = None
        nc = 0.0
        try:
            nc = float(fObj.coefficient)
        except ValueError:
            print 'Suspected rational number (%s) detected in fluxObjective %s' % (fObj.coefficient, fObj.getPid())
        if nc >= 0.0:
            sign = '+'
        else:
            sign = '-'    
        # TODO: if use_rational is not used simply try and evaluate the coefficient string with
        # sympy.Rational.evalf()  and use this as the value for nc ... also remove use_rational case
        if use_rational:
            objStr += ' %s%s%s' % (sympy.Rational(fObj.coefficient), multisymb, fObj.reaction)
            ##  FF.write('%s: %s%s%s\n' % (fobj0.reaction, sympy.Rational(fobj0.coefficient), multisymb, fobj0.reaction))
        else:
            objStr += ' %s %s%s%s' % (sign, abs(nc), multisymb, fObj.reaction)
            ##  FF.write('%s: %s%s%s\n' % (fobj0.reaction, fobj0.coefficient, multisymb, fobj0.reaction))
    FF.write('%s\n' % objStr)
    if constraint_mode == 'math':
        CONST = BuildLPConstraintsMath(fba)
    elif constraint_mode == 'relaxed':
        CONST = BuildLPConstraintsRelaxed(fba)
    else:
        CONST = BuildLPConstraintsStrict(fba, use_rational)

    CONST.seek(0)
    BOUNDS = BuildLPFluxBounds(fba, use_rational)
    BOUNDS.seek(0)
    if __DEBUG__: print CONST.read(); CONST.seek(0)
    if __DEBUG__: print BOUNDS.read(); BOUNDS.seek(0)

    FF.write('\nSubject To\n')
    FF.write(CONST.read())
    FF.write('\nBounds\n')
    FF.write(BOUNDS.read())
    FF.write('\nEND\n')
    FF.close()
    print 'WriteModelLP has written a file to %s' % FNAME
    return FNAME

def WriteModelLP(fba, work_dir=None, fname=None, multisymb=' ', format='%s', use_rational=False, constraint_mode=None, quiet=False):
    """
    Writes an FBA object as an LP in CPLEX LP format

     - *fba* an instantiated FBAmodel instance
     - *work_dir* directory designated for output
     - *fname* the file name [default=fba.getPid()]
     - *multisymb* the multiplication symbol (default: <space>)
     - *format* the number format of the output
     - *use_rational* output rational numbers [default=False]
     - *quiet* [default=False] supress information messages

    """
    if constraint_mode != None:
        print "\nConstraint_mode has been deprecated"
    
    if not _HAVE_SYMPY_ and use_rational:
        use_rational = False
        print '\nWarning switching to floating point arithmetic: install Sympy for rational IO'    
        time.sleep(2)
    
    FNAME = None
    if fname == None:
        fname = fba.getPid()
    if work_dir != None:
        FnameTmp = os.path.join(work_dir,  fname)
    else:
        FnameTmp = fname

    FNAME = FnameTmp+'.lp'
    FF = file(FNAME, 'w')
    FF.write('\\\\ %s \n\n' % FnameTmp)
    if len(fba.objectives) > 0:
        ##  print fba.objectives
        if fba.objectives[fba.activeObjIdx].operation == None:
            print '\nWARNING: Objective function \"%s\" has no \"operation\" defined assuming \"maximize\"' % fba.objectives[fba.activeObjIdx].getPid()
            fba.objectives[fba.activeObjIdx].operation = 'maximize'
            time.sleep(2)
        objO = fba.objectives[fba.activeObjIdx].operation.lower()
        objO = objO[0].upper() + objO[1:]
        FF.write('%s\n' % objO)
        objStr = '%s_objf: ' % fba.objectives[fba.activeObjIdx].getPid()
        for fObj in fba.objectives[fba.activeObjIdx].fluxObjectives:
            sign = None
            nc = 0.0
            try:
                nc = float(fObj.coefficient)
            except ValueError:
                if _HAVE_SYMPY_:
                    nc = sympy.Rational(fObj.coefficient).evalf()
                else:
                    raise ValueError, 'Invalid coefficient (%s) detected in fluxObjective %s' % (fObj.coefficient, fObj.getPid())
            if nc >= 0.0:
                sign = '+'
            else:
                sign = '-'
            if use_rational:
                objStr += ' %s %s%s%s' % (sign, sympy.Rational(abs(nc)), multisymb, fObj.reaction)
            else:
                objStr += ' %s %s%s%s' % (sign, abs(nc), multisymb, fObj.reaction)
    else:
        objStr = '\n\\\\ No objectives defined\n\n'

    FF.write('%s\n' % objStr)
    CONST = BuildLPConstraints(fba, use_rational)
    CONST.seek(0)
    if fba.user_constraints != None and len(fba.user_constraints.keys()) > 0:
        UCONST = BuildLPUserConstraints(fba, use_rational=False)
        UCONST.seek(0)
    BOUNDS = BuildLPFluxBounds(fba, use_rational)
    BOUNDS.seek(0)
    if __DEBUG__: print CONST.read(); CONST.seek(0)
    if __DEBUG__: print BOUNDS.read(); BOUNDS.seek(0)

    FF.write('\nSubject To\n')
    FF.write(CONST.read())
    if fba.user_constraints != None and len(fba.user_constraints.keys()) > 0:
        FF.write('\\\\UserConstraints\n')
        FF.write(UCONST.read())
    FF.write('\nBounds\n')
    FF.write(BOUNDS.read())
    FF.write('\nEND\n')
    FF.close()
    if not quiet:
        print 'WriteModelLP has written a file to %s' % FNAME
    return FNAME

def BuildHformatFluxBounds(fba, infinity_replace=None):
    """
    Build and return a cStringIO that contains the flux bounds in H format
    
     - *fba* a PySCeS-CBM FBA object
     - *infinity_replace* [default=None] if defined this is the abs(value) of +-<infinity>
    
    """
    LBs = {}
    UBs = {}
    for c in fba.flux_bounds:
        ##  print '%s: %s' % (c.getPid(), c.operation)
        minv = None
        maxv = None
        oper = None
        R = c.reaction
        if c.operation in ['less','lessEqual']:
            if infinity_replace != None and numpy.isposinf([c.value])[0]:
                maxv = infinity_replace
            elif infinity_replace != None and numpy.isneginf([c.value])[0]:
                maxv = -infinity_replace
            else:
                maxv = c.value
            ##  print 'maxv', maxv                
            oper = c.operation
            if oper == 'less':
                oper = '<='
                ##  print 'LP Bounds operator must be <= not <'
            else:
                oper = '<='
        if c.operation in ['greater','greaterEqual']:
            if infinity_replace != None and numpy.isposinf([c.value])[0]:
                minv = infinity_replace
            elif infinity_replace != None and numpy.isneginf([c.value])[0]:
                minv = -infinity_replace
            else:
                minv = c.value
            ##  print 'minv', minv                
            oper = c.operation
            if oper == 'greater':
                oper = '<='
                ##  print 'LP Bounds operator must be <= not <'
            else:
                oper = '<='
        if c.operation == 'equals':
            if infinity_replace != None and numpy.isposinf([c.value])[0]:
                maxv = infinity_replace
            elif infinity_replace != None and numpy.isneginf([c.value])[0]:
                maxv = -infinity_replace
            else:
                maxv = c.value                        
            oper = '='
            ##  print 'maxv2', maxv                
        if maxv != None:
            UBs.update({R : float(maxv)})
        elif minv != None:
            LBs.update({R : float(minv)})

    if __DEBUG__:
        print ' '
        print LBs
        print UBs

    BsRHS = []
    LBm = numpy.zeros((len(LBs), fba.N.shape[1]))
    UBm = numpy.zeros((len(UBs), fba.N.shape[1]))

    for lb in range(len(LBs.keys())):
        LBm[lb, fba.N.col.index(LBs.keys()[lb])] = 1.0
        BsRHS.append(LBs[LBs.keys()[lb]])

    if __DEBUG__:
        print fba.N.col
        print LBm
        print UBm
        print BsRHS

    for ub in range(len(UBs.keys())):
        UBm[ub, fba.N.col.index(UBs.keys()[ub])] = -1.0
        BsRHS.append(-UBs[UBs.keys()[ub]])

    if __DEBUG__:
        print fba.N.col
        print LBm
        print UBm
        print BsRHS
        print ' '

    return numpy.vstack([LBm, UBm]), BsRHS

def WriteModelHFormatFBA(fba, work_dir=None, use_rational=False, fullLP=True, format='%s', infinity_replace=None):
    """
    Write an FBA-LP in polynomial H-Format file. This version has been replaced byt `WriteModelHFormatFBA2()`
    but is kept for backwards compatability.
    
     - *fba* a PySCeS-CBM FBA object
     - *Work_dir* [default=None] the output directory
     - *use_rational* [default=false] use rational numbers in output (requires sympy)
     - *fullLP* [default=True] include the default objective function as a maximization target
     - *format* [default='%s'] the number format string
     - *infinity_replace* [default=None] if defined this is the abs(value) of +-<infinity> 
    
    """
    
    print '\nThis method is deprecated please use: WriteModelHFormatFBA2\n'
    
    if not _HAVE_SYMPY_ and use_rational:
        use_rational = False
        print 'Warning: install Sympy for rational IO'
    M = fba
    LHS = M.N.array.copy()
    RHS = [0.0 for e in range(M.N.shape[0])]
    if __DEBUG__:
        print LHS
        print RHS
    LHS = numpy.vstack([LHS, -M.N.array.copy()])
    RHS += [0.0 for e in range(M.N.shape[0])]
    if __DEBUG__:
        print LHS
        print RHS
    BsLHS, BsRHS = BuildHformatFluxBounds(M, infinity_replace=infinity_replace)
    if __DEBUG__:
        print BsLHS
        print BsRHS
    LHS = numpy.vstack([LHS, BsLHS])
    RHS += BsRHS
    del BsLHS, BsRHS
    if __DEBUG__:
        print LHS
        print RHS
    if not use_rational:
        name = M.getPid().replace('.xml', '') + '.ine'
    else:
        name = M.getPid().replace('.xml', '') + '_r.ine'

    RHS = numpy.array(RHS,'d')
    RHS.shape = (len(RHS), 1)
    if __DEBUG__: print RHS
    ##  LP = numpy.hstack([LHS, RHS])

    OBJ_FUNC = numpy.zeros(LHS.shape[1]+1)
    for j in range(LHS.shape[1]):
        # first objective function, first flux objective
        if __DEBUG__: print M.objectives[0].fluxObjectives[0].reaction, M.N.col[j]
        if M.objectives[0].fluxObjectives[0].reaction == M.N.col[j]:
            OBJ_FUNC[j] = float(M.objectives[0].fluxObjectives[0].coefficient)
    if __DEBUG__: print OBJ_FUNC

    # for Ax >= B Hformat wants -B A >= 0
    LP = numpy.hstack([-RHS, LHS])
    OBJ_FUNC = numpy.hstack([-OBJ_FUNC[-1], OBJ_FUNC[:-1]])

    if __DEBUG__:
        print OBJ_FUNC
        print LP
    del LHS, RHS

    if work_dir == None:
        Fname = name
    else:
        assert os.path.exists(work_dir), '\nJanee ...'
        Fname = os.path.join(work_dir, name)

    F = file(Fname, 'w')
    F.write('* %s\n\n' % name)
    F.write('H-representation\n\nbegin\n')
    NUM_TYPE = 'real'
    if use_rational:
        NUM_TYPE = 'rational'
    F.write('%s  %s  %s\n' % (LP.shape[0], LP.shape[1], NUM_TYPE))

    strW = format+' '
    for r in range(LP.shape[0]):
        for c in range(LP.shape[1]):
            if not use_rational:
                if LP[r,c] == 0.0 or LP[r,c] == -0.0:
                    LP[r,c] = 0.0
                F.write(strW % LP[r,c])
            else:
                F.write('%s ' % sympy.Rational(format % LP[r,c]))
        F.write('\n')

    if fullLP:
        F.write('end\nlponly\n')
        F.write('maximize\n')
        for o in OBJ_FUNC:
            if not use_rational:
                F.write(strW % o)
            else:
                F.write('%s ' % sympy.Rational(format % o))
    else:
        F.write('end\n')
    F.write('\n')
    F.close()
    F = file(Fname.replace('.ine','')+'.columns.txt', 'w')
    for j in range(M.N.array.shape[1]):
        F.write('%s,%s\n' % (j, M.N.col[j]))
    F.write('\n')
    F.close()
    return Fname

def WriteModelHFormatFBA2(fba, fname=None, work_dir=None, use_rational=False, fullLP=True, format='%s', infinity_replace=None):
    """
    Write an FBA-LP in polynomial H-Format file. This is an improved version of `WriteModelHFormatFBA()`
    which it replaces but is kept for backwards compatability.
    
     - *fba* a PySCeS-CBM FBA object
     - *fname* [default=None] the output filename, fba.getPid() if not defined
     - *Work_dir* [default=None] the output directory
     - *use_rational* [default=false] use rational numbers in output (requires sympy)
     - *fullLP* [default=True] include the default objective function as a maximization target
     - *format* [default='%s'] the number format string
     - *infinity_replace* [default=None] if defined this is the abs(value) of +-<infinity>
    
    """
    if not _HAVE_SYMPY_ and use_rational:
        use_rational = False
        print 'Warning: install Sympy for rational IO'
    M = fba
    LHS = M.N.array.copy()
    RHS = [0.0 for e in range(M.N.shape[0])]
    if __DEBUG__:
        print LHS
        print RHS
    LHS = numpy.vstack([LHS, -M.N.array.copy()])
    RHS += [0.0 for e in range(M.N.shape[0])]
    if __DEBUG__:
        print LHS
        print RHS
    BsLHS, BsRHS = BuildHformatFluxBounds(M, infinity_replace=infinity_replace)
    if __DEBUG__:
        print BsLHS
        print BsRHS
    LHS = numpy.vstack([LHS, BsLHS])
    RHS += BsRHS
    del BsLHS, BsRHS
    if __DEBUG__:
        print LHS
        print RHS

    RHS = numpy.array(RHS,'d')
    RHS.shape = (len(RHS), 1)
    if __DEBUG__: print RHS
    ##  LP = numpy.hstack([LHS, RHS])

    OBJ_FUNC = numpy.zeros(LHS.shape[1]+1)
    objIdx = M.activeObjIdx
    for j in range(LHS.shape[1]):
        for fo in range(len(M.objectives[objIdx].getFluxObjectiveReactions())):
            if M.objectives[objIdx].fluxObjectives[fo].reaction == M.N.col[j]:
                print M.objectives[objIdx].fluxObjectives[fo].reaction, M.N.col[j]
                OBJ_FUNC[j] = float(M.objectives[objIdx].fluxObjectives[fo].coefficient)
    ##  print OBJ_FUNC

    # for Ax >= B Hformat wants -B A >= 0
    LP = numpy.hstack([-RHS, LHS])
    OBJ_FUNC = numpy.hstack([-OBJ_FUNC[-1], OBJ_FUNC[:-1]])

    if __DEBUG__:
        print OBJ_FUNC
        print LP
    del LHS, RHS

    if work_dir != None:
        assert os.path.exists(work_dir), '\nJanee ...'
        fname = os.path.join(work_dir, fname)

    if fname == None:
        fname = M.getPid().replace('.xml', '')
    if not use_rational:
        fname += '.ine'
    else:
        fname += '_r.ine'

    F = file(fname, 'w')
    F.write('* %s\n' % os.path.split(fname)[-1])
    F.write('H-representation\nbegin\n')
    NUM_TYPE = 'real'
    if use_rational:
        NUM_TYPE = 'rational'
    F.write('%s  %s  %s\n' % (LP.shape[0], LP.shape[1], NUM_TYPE))

    strW = format+' '
    for r in range(LP.shape[0]):
        for c in range(LP.shape[1]):
            if not use_rational:
                if LP[r,c] == 0.0 or LP[r,c] == -0.0:
                    LP[r,c] = 0.0
                F.write(strW % LP[r,c])
            else:
                ##  print LP[r,c]
                F.write('%s ' % sympy.Rational(format % LP[r,c]))
        F.write('\n')

    if fullLP:
        F.write('end\nlponly\n')
        F.write('maximize\n') # check if Hformat has a minimize kw
        for o in OBJ_FUNC:
            if not use_rational:
                F.write(strW % o)
            else:
                F.write('%s ' % sympy.Rational(format % o))        
        # then we can use use this
        ##  F.write('%s\n' % M.objectives[M.activeObjIdx].operation)
        ##  if M.activeObjIdx].operation == 'maximize':
            ##  for o in OBJ_FUNC:
                ##  if not use_rational:
                    ##  F.write(strW % o)
                ##  else:
                    ##  F.write('%s ' % sympy.Rational(format % o))
        ##  else:
            ##  for o in OBJ_FUNC:
                ##  o = -o
                ##  if not use_rational:
                    ##  F.write(strW % o)
                ##  else:
                    ##  F.write('%s ' % sympy.Rational(format % o))            
    else:
        F.write('end\n')
    F.write('\n')
    F.close()
    F = file(fname.replace('.ine','')+'.columns.txt', 'w')
    for j in range(M.N.array.shape[1]):
        F.write('%s,%s\n' % (j, M.N.col[j]))
    F.write('\n')
    F.close()
    return fname

def writeListToLP(fname, obj=None, const=None, bnds=None, work_dir=None, objtype='maximize'):
    if work_dir == None:
        work_dir = os.getcwd()
    F = file(os.path.join(work_dir, fname+'.lp'), 'w')
    F.write("\\\\ %s\n" % fname)
    objtype = objtype.lower()
    if objtype == 'max': objtype = 'maximize'
    if objtype == 'min': objtype = 'minimize'    
    if objtype in ['maximise', 'minimise']:
        objtype = objtype.replace('se','ze')
    assert objtype in ['maximize', 'minimize'], "\nobjtype must be ['maximize', 'minimize'] not %s" % objtype  

    if obj != None:
        if objtype == 'maximize':
            F.write('\nMaximize\n')
        elif objtype == 'minimize':
            F.write('\nMinimize\n')
        for o in obj:
            F.write(' %s\n' % o)
    if const != None:
        F.write('\nSubject to\n')
        for c in const:
            F.write('%s\n' % c)
    if bnds != None:
        F.write('\nBounds\n')
        for b in bnds:
            F.write('%s\n' % b)
    F.write('\nEND\n')
    F.close()
    print 'LP written to: %s.lp' % os.path.join(work_dir, fname)
    return os.path.join(work_dir, fname+'.lp')


def writeMinDistanceLP(fname, fbas, work_dir=None, ignoreDistance=[], with_protein_cost=False, constraint_mode='strict'):
    if work_dir == None:
        work_dir = os.getcwd()

    fC = []
    objFname = ''
    for l in fbas:
        fC.append(len(l.reactions))
        objFname += l.prefix
    fC = numpy.array(fC)
    if __DEBUG__:
        print fC
        print (fC == fC[0])
        print numpy.alltrue((fC == fC[0]))
    assert numpy.alltrue((fC == fC[0])), '\nModels must have the same number of fluxes\n!'

    conL = []
    # model flux_bounds
    initial_cnstr = []

    for f in fbas:
        if constraint_mode == 'math':
            initial_cnstr.append(BuildLPConstraintsMath(f))
        elif constraint_mode == 'relaxed':
            initial_cnstr.append(BuildLPConstraintsRelaxed(f))
        else:
            initial_cnstr.append(BuildLPConstraintsStrict(f))

    for ib in initial_cnstr:
        ib.seek(0)
        for l in ib:
            conL.append(l.strip())
        conL.append(' ')
    del initial_cnstr

    bndL = []
    # model bounds
    initial_bnds = []
    for f in fbas:
        initial_bnds.append(BuildLPFluxBounds(f))
    for ib in initial_bnds:
        ib.seek(0)
        for l in ib:
            bndL.append(l.strip())
        bndL.append(' ')
    del initial_bnds

    artVar = []
    artVarX = []
    ignoreDistance = []

    Combi = CBTools.ComboGen()
    Cnumber = 2
    Cdata = ''
    unique_combinations = None
    for x in range(len(fbas)):
        Cdata += '%s' % x
    if __DEBUG__: print Cdata

    Combi.uniqueCombinations(Cdata, Cnumber, temp=[])
    Combi.numberifyComb2Int()
    unique_combinations = Combi.combo_int

    if __DEBUG__:
        print 'Data (%s):\n%s\n' % (Cnumber, Cdata)
        print 'UniqueCombinStr:\n%s' % Combi.combo
        print 'UniqueCombinations:\n%s' % unique_combinations

    zbase = 0
    combcount = 1
    for uq in unique_combinations:
        if __DEBUG__: print uq
        for s in range(len(fbas[uq[0]].reactions)):
            RiD1 = fbas[uq[0]].getReactionIds()
            RiD2 = fbas[uq[1]].getReactionIds()
            if __DEBUG__: print RiD1[s], RiD2[s]
            if RiD1[s] not in ignoreDistance:
                ##  av = 'z%s' % (zbase+s+1)
                av = 'zvar%s%s' % (combcount,RiD1[s].replace(fbas[uq[0]].prefix,''))
                
                c1 = '%sa: %s - %s - %s <= 0.0' % (av, RiD1[s], RiD2[s], av)
                ##  c1 = '%s - %s - %s <= 0.0' % (RiD1[s], RiD2[s], av)
                c2 = '%sb: %s - %s + %s >= 0.0' % (av, RiD1[s], RiD2[s], av)
                ##  c2 = '%s - %s + %s >= 0.0' % (RiD1[s], RiD2[s], av)
                
                # add the protein cost 
                if with_protein_cost:
                    av = '%s %s' % (fbas[uq[0]].reactions[s].annotation['CBM_PEPTIDE_COST'], av)
                artVar.append(av)
                conL.append(c1)
                conL.append(c2)

        combcount += 1
        zbase += len(fbas[uq[0]].reactions)

    objS = '%smulti: ' % objFname
    vcntr = 0
    for v in artVar:
        objS += '%s + ' % v
        vcntr += 1
        if vcntr >= 500:
            objS += '\n'
            vcntr = 0
    objS = objS[:-3]
    if len(artVarX) >= 1:
        objS += ' \\* Ignored: '
        for o in artVarX:
            objS += '%s ' % o
        objS += '*\\\n'
    objL = [objS]

    objFcnstr = [' ']
    assert len(f.objectives[f.activeObjIdx].getFluxObjectiveReactions()) == 1, "\nOnly single fluxObjectives dealt with at this time"
    for f in fbas:
        objFcnstr.append('C_%s: %s >= %f' % (f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0],\
                                    f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0], f.objectives[f.activeObjIdx].value))
    conL = conL + objFcnstr

    if __DEBUG__:
        print objL
        for c in conL:
            print c
        print ' '
        for b in bndL:
            print b

    F = file(os.path.join(work_dir, fname+'.lp'), 'w')
    header = '\\\\ MultiInputMinimization: '
    for f in fbas:
        header += '%s, ' % f.getPid()
    F.write('%s\n' % header[:-2])
    F.write('\nMinimize\n')
    for o in objL:
        F.write('%s \n' % o)
    F.write('\nSubject to\n')
    for c in conL:
        F.write('%s \n' % c)
    F.write('\nBounds\n')
    for b in bndL:
        F.write(' %s \n' % b)
    F.write('END\n\n')
    F.close()
    print 'LP written to: %s.lp' % os.path.join(work_dir, fname)
    return os.path.join(work_dir, fname+'.lp')

def writeMinDistanceLP_absL1(fname, fbas, work_dir=None, ignoreDistance=[], bigM=500, with_protein_cost=False, constraint_mode='strict'):
    if work_dir == None:
        work_dir = os.getcwd()

    fC = []
    objFname = ''
    for l in fbas:
        fC.append(len(l.reactions))
        objFname += l.prefix
    fC = numpy.array(fC)
    if __DEBUG__:
        print fC
        print (fC == fC[0])
        print numpy.alltrue((fC == fC[0]))
    assert numpy.alltrue((fC == fC[0])), '\nModels must have the same number of fluxes\n!'

    conL = []
    # model flux_bounds
    initial_cnstr = []

    for f in fbas:
        if constraint_mode == 'math':
            initial_cnstr.append(BuildLPConstraintsMath(f))
        elif constraint_mode == 'relaxed':
            initial_cnstr.append(BuildLPConstraintsRelaxed(f))
        else:
            initial_cnstr.append(BuildLPConstraintsStrict(f))

    for ib in initial_cnstr:
        ib.seek(0)
        for l in ib:
            conL.append(l.strip())
        conL.append(' ')
    del initial_cnstr

    bndL = []
    # model bounds
    initial_bnds = []
    bigMS = []
    for f in fbas:
        initial_bnds.append(BuildLPFluxBounds(f))
        ##  bigMS.append(max([abs(float(v.value)) for v in f.reactions]))
    for ib in initial_bnds:
        ib.seek(0)
        for l in ib:
            bndL.append(l.strip())
        bndL.append(' ')
    del initial_bnds

    Combi = CBTools.ComboGen()
    Cnumber = 2
    Cdata = ''
    unique_combinations = None
    for x in range(len(fbas)):
        Cdata += '%s' % x
    if __DEBUG__: print Cdata

    Combi.uniqueCombinations(Cdata, Cnumber, temp=[])
    Combi.numberifyComb2Int()
    unique_combinations = Combi.combo_int

    if __DEBUG__:
        print 'Data (%s):\n%s\n' % (Cnumber, Cdata)
        print 'UniqueCombinStr:\n%s' % Combi.combo
        print 'UniqueCombinations:\n%s' % unique_combinations
    print bigMS
    print '\nbigM = ', bigM, '\n'
    
    
    ##  ILPMETHOD = 'SK' # steven
    ILPMETHOD = 'GK' # gunnar    
    zbase = 0
    combcount = 1
    
    artVar = []
    artVarX = []
    ignoreDistance = []    
    boolVars = []

    for uq in unique_combinations:
        if __DEBUG__: print uq
        for s in range(len(fbas[uq[0]].reactions)):
            RiD1 = fbas[uq[0]].getReactionIds()
            RiD2 = fbas[uq[1]].getReactionIds()
            if __DEBUG__: print RiD1[s], RiD2[s]
            if RiD1[s] not in ignoreDistance:
                Var1 = RiD1[s]
                bVar1 =  'xvar_%s' % Var1
                absVar1 = 'absL_%s' % Var1
                Var2 = RiD2[s]
                bVar2 =  'xvar_%s' % Var2
                absVar2 = 'absL_%s' % Var2
                
                c0a = '\n'                
                c0a += '%s - %s >= 0\n' % (absVar1, Var1)
                c0a += '%s + %s >= 0\n' % (absVar1, Var1)
                
                if ILPMETHOD == 'GK':
                    # gunnar
                    c0a += '%s + %s - %s %s <= 0\n' % (absVar1, Var1, bigM, bVar1)
                    c0a += '%s - %s + %s %s <= %s\n' % (absVar1, Var1, bigM, bVar1, bigM)
                elif ILPMETHOD == 'SK':
                    # steven
                    c0a += '%s - %s - %s %s <= 0\n' % (absVar1, Var1, bigM, bVar1)
                    c0a += '%s + %s + %s %s <= %s\n' % (absVar1, Var1, bigM, bVar1, bigM)
                    
                c0b = '\n'                
                c0b += '%s - %s >= 0\n' % (absVar2, Var2)
                c0b += '%s + %s >= 0\n' % (absVar2, Var2)
                if ILPMETHOD == 'GK':
                    # gunnar
                    c0b += '%s + %s - %s %s <= 0\n' % (absVar2, Var2, bigM, bVar2)
                    c0b += '%s - %s + %s %s <= %s\n' % (absVar2, Var2, bigM, bVar2, bigM)
                elif ILPMETHOD == 'SK':
                    # steven
                    c0b += '%s - %s - %s %s <= 0\n' % (absVar2, Var2, bigM, bVar2)
                    c0b += '%s + %s + %s %s <= %s\n' % (absVar2, Var2, bigM, bVar2, bigM)

                c0 = c0a + c0b

                av = 'zvar%s%s' % (combcount,Var1.replace(fbas[uq[0]].prefix,''))
                
                if bVar1 not in boolVars:
                    boolVars.append(bVar1)
                if bVar2 not in boolVars:
                    boolVars.append(bVar2)
                c1 = '%sa: %s - %s - %s <= 0.0' % (av, absVar1, absVar2, av)
                c2 = '%sb: %s - %s + %s >= 0.0' % (av, absVar1, absVar2, av)
                
                # add the protein cost 
                if with_protein_cost:
                    av = '%s %s' % (fbas[uq[0]].reactions[s].annotation['CBM_PEPTIDE_COST'], av)
                artVar.append(av)
                
                conL.append(c0)
                conL.append(c1)
                conL.append(c2)
        combcount += 1
        zbase += len(fbas[uq[0]].reactions)

    objS = '%smulti: ' % objFname
    vcntr = 0
    for v in artVar:
        objS += '%s + ' % v
        vcntr += 1
        if vcntr >= 500:
            objS += '\n'
            vcntr = 0
    objS = objS[:-3]
    if len(artVarX) >= 1:
        objS += ' \\* Ignored: '
        for o in artVarX:
            objS += '%s ' % o
        objS += '*\\\n'
    objL = [objS]

    objFcnstr = [' ']
    assert len(f.objectives[f.activeObjIdx].getFluxObjectiveReactions()) == 1, "\nOnly single fluxObjectives dealt with at this time"
    for f in fbas:
        OFvalue = f.objectives[f.activeObjIdx].value
        objFcnstr.append('C_%s: %s >= %s' % (f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0],\
                                    f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0], OFvalue))
    conL = conL + objFcnstr

    if __DEBUG__:
        print objL
        for c in conL:
            print c
        print ' '
        for b in bndL:
            print b

    F = file(os.path.join(work_dir, fname+'.lp'), 'w')
    header = '\\\\ MultiInputMinimization: '
    for f in fbas:
        header += '%s, ' % f.getPid()
    F.write('%s\n' % header[:-2])
    F.write('\nMinimize\n')
    for o in objL:
        F.write('%s \n' % o)
    F.write('\nSubject to\n')
    for c in conL:
        F.write('%s \n' % c)
    F.write('\nBounds\n')
    for b in bndL:
        F.write(' %s \n' % b)
    if len(boolVars) > 0:
        F.write('Binary\n')
        for b in boolVars:
            F.write(' %s \n' % b)        
    F.write('\nEND\n')
    F.close()
    print 'LP written to: %s.lp' % os.path.join(work_dir, fname)
    return os.path.join(work_dir, fname+'.lp')



def writeMinDistanceLPwithCost(fname, fbas, work_dir=None, ignoreDistance=[], constraint_mode='strict'):
    """
    For backwards compatability only
    """
    print "\n\n**********\nDeprecation warning!\nPlease use writeMinDistanceLP(with_protein_cost=True) instead of writeMinDistanceLPwithCost()\n\n**********\n"
    writeMinDistanceLP(fname, fbas, work_dir=work_dir, ignoreDistance=ignoreDistance, with_protein_cost=True, constraint_mode=constraint_mode)

# CAN GO SOON JUST PUTTING IN FOR SVN SYNCH
"""
def writeMinDistanceLPwithCost(fname, fbas, work_dir=None, ignoreDistance=[], constraint_mode='strict'):
    if work_dir == None:
        work_dir = os.getcwd()

    fC = []
    objFname = ''
    for l in fbas:
        fC.append(len(l.reactions))
        objFname += l.prefix
    fC = numpy.array(fC)
    if __DEBUG__:
        print fC
        print (fC == fC[0])
        print numpy.alltrue((fC == fC[0]))
    assert numpy.alltrue((fC == fC[0])), '\nModels must have the same number of fluxes\n!'

    conL = []
    # model flux_bounds
    initial_cnstr = []

    for f in fbas:
        if constraint_mode == 'math':
            initial_cnstr.append(BuildLPConstraintsMath(f))
        elif constraint_mode == 'relaxed':
            initial_cnstr.append(BuildLPConstraintsRelaxed(f))
        else:
            initial_cnstr.append(BuildLPConstraintsStrict(f))

    for ib in initial_cnstr:
        ib.seek(0)
        for l in ib:
            conL.append(l.strip())
        conL.append(' ')
    del initial_cnstr

    bndL = []
    # model bounds
    initial_bnds = []
    for f in fbas:
        initial_bnds.append(BuildLPFluxBounds(f))
    for ib in initial_bnds:
        ib.seek(0)
        for l in ib:
            bndL.append(l.strip())
        bndL.append(' ')
    del initial_bnds

    artVar = []
    artVarX = []
    ignoreDistance = []

    Combi = CBTools.ComboGen()
    Cnumber = 2
    Cdata = ''
    unique_combinations = None
    for x in range(len(fbas)):
        Cdata += '%s' % x
    if __DEBUG__: print Cdata

    Combi.uniqueCombinations(Cdata, Cnumber, temp=[])
    Combi.numberifyComb2Int()
    unique_combinations = Combi.combo_int

    if __DEBUG__:
        print 'Data (%s):\n%s\n' % (Cnumber, Cdata)
        print 'UniqueCombinStr:\n%s' % Combi.combo
        print 'UniqueCombinations:\n%s' % unique_combinations

    zbase = 0
    combcount = 1
    for uq in unique_combinations:
        if __DEBUG__: print uq
        RiD1 = fbas[uq[0]].getReactionIds()
        RiD2 = fbas[uq[1]].getReactionIds()
        for s in range(len(fbas[uq[0]].reactions)):
            if __DEBUG__: print RiD1[s], RiD2[s]
            if RiD1[s] not in ignoreDistance:
                ##  av = 'z%s' % (zbase+s+1)
                ##  av = 'zvar%s' % (RiD1[s].replace(fbas[uq[0]].prefix,''))
                av = 'zvar%s%s' % (combcount,RiD1[s].replace(fbas[uq[0]].prefix,''))
                c1 = '%sa: %s - %s - %s <= 0.0' % (av, RiD1[s], RiD2[s], av)
                ##  c1 = '%s - %s - %s <= 0.0' % (RiD1[s], RiD2[s], av)
                c2 = '%sb: %s - %s + %s >= 0.0' % (av, RiD1[s], RiD2[s], av)
                ##  c2 = '%s - %s + %s >= 0.0' % (RiD1[s], RiD2[s], av)
                conL.append(c1)
                conL.append(c2)
                ##  print RiD1[s], RiD2[s], av
                ##  print fbas[uq[0]].reactions[s].getPid(), fbas[uq[1]].reactions[s].getPid()
                ##  print fbas[uq[0]].reactions[s].annotation['CBM_PEPTIDE_COST'], fbas[uq[1]].reactions[s].annotation['CBM_PEPTIDE_COST']
                av = '%s %s' % (fbas[uq[0]].reactions[s].annotation['CBM_PEPTIDE_COST'], av)
                ##  print RiD1[s], RiD2[s], av
                artVar.append(av)
            else:
                # THIS MAY BE A BUG
                ##  av = 'z%s' % (zbase+s+1)
                ##  av = 'zvar%s' % (RiD1[s].replace(fbas[uq[0]].prefix,''))
                av = 'zvar%s%s' % (combcount,RiD1[s].replace(fbas[uq[0]].prefix,''))
                artVarX.append(av)
                c1 = '\\* %sa: %s - %s - %s <= 0.0 *\\' % (av, RiD1[s], RiD2[s], av)
                ##  c1 = '\\* %s - %s - %s <= 0.0 *\\' % (RiD1[s], RiD2[s], av)
                c2 = '\\* %sb: %s - %s + %s >= 0.0 *\\' % (av, RiD1[s], RiD2[s], av)
                ##  c2 = '\\* %s - %s + %s >= 0.0 *\\' % (RiD1[s], RiD2[s], av)
                conL.append(c1)
                conL.append(c2)
        combcount += 1
        zbase += len(fbas[uq[0]].reactions)

    objS = '%smulti: ' % objFname
    vcntr = 0
    for v in artVar:
        objS += '%s + ' % v
        vcntr += 1
        if vcntr >= 500:
            objS += '\n'
            vcntr = 0
    objS = objS[:-3]
    if len(artVarX) >= 1:
        objS += ' \\* Ignored: '
        for o in artVarX:
            objS += '%s ' % o
        objS += '*\\\n'
    objL = [objS]

    objFcnstr = [' ']
    for f in fbas:
        objFcnstr.append('C_%s: %s >= %f' % (f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0],\
                                    f.objectives[f.activeObjIdx].getFluxObjectiveReactions()[0], f.objectives[f.activeObjIdx].value))
    conL = conL + objFcnstr

    if __DEBUG__:
        print objL
        for c in conL:
            print c
        print ' '
        for b in bndL:
            print b

    F = file(os.path.join(work_dir, fname+'.lp'), 'w')
    header = '\\\\ MultiInputMinimization: '
    for f in fbas:
        header += '%s, ' % f.getPid()
    F.write('%s\n' % header[:-2])
    F.write('\nMinimize\n')
    for o in objL:
        F.write('%s \n' % o)
    F.write('\nSubject to\n')
    for c in conL:
        F.write('%s \n' % c)
    F.write('\nBounds\n')
    for b in bndL:
        F.write(' %s \n' % b)
    F.write('END\n\n')
    F.close()
    print 'LP written to: %s.lp' % os.path.join(work_dir, fname)
    return os.path.join(work_dir, fname+'.lp')
"""

def writeOptimalSolution(fba, fname, Dir=None, separator=',', only_exchange=False):
    """
    This function writes the optimal solution to file
    
     - *fba* an instance of an PySCeSCBM model
     - *fname* the output filename
     - *Dir* [default=None] use directory if not None
     - *separator* [default=','] the column separator
     - *only_exchange* [default=False] only output fluxes labelled as exchange reactions
     
    """    
    if Dir != None:
        assert os.path.exists(Dir), '\nPath does not exist'
        fname = os.path.join(Dir, fname)
    fname_r = fname + '_solution.csv'
    objName = ''
    if len(fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions()) > 1:
        for J in fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions():
            objName += '%s_' % J
        objName = objName[:-1]
    else:
        objName = fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions()[0]
    try:
        F = file(fname_r, 'w')
    except IOError:
        print '\nCSV file \"%s\" is locked by an external application (probably Excel) please close file and try again (or use a different filename).' % fname_r
        return 
    cntr = 0
    F.write('%s%s%s%s%s%s\n' % ('ObjectiveFunction',separator,objName,separator,separator,separator))
    F.write('\"%s\"%s%s%s%s%s%s%s%s%s\"%s\"%s\"%s\"\n' % ('Reaction',separator,'Value',separator,'LowerBound',separator,'UpperBound',separator,'Reduced cost',separator,'Name',separator,'Gene association'))
    for r in fba.reactions:
        GO = False
        if not only_exchange:
            GO = True
        elif only_exchange and r.is_exchange:
            GO = True
        if GO:
            if r.annotation.has_key('GENE ASSOCIATION'):
                gene = r.annotation['GENE ASSOCIATION']
            else:
                gene = 'none'
            bnds = fba.getReactionBounds(r.getPid())
            Lbnd = -numpy.inf
            Ubnd = numpy.inf
            if bnds != None:
                if bnds[1] != None:
                    Lbnd = bnds[1]
                if bnds[2] != None:
                    Ubnd = bnds[2]
            F.write('\"%s\"%s%s%s%s%s%s%s%s%s\"%s\"%s\"%s\"\n' % (r.getPid(),separator,r.value,separator,Lbnd,separator,Ubnd,separator,r.reduced_cost,separator,r.name,separator,gene))
    F.flush()
    F.close()
    print 'Reactions exported as CSV to %s' % fname_r

def writeModelToCSV(fba, fname, Dir=None, separator=',', only_exchange=False):
    """
    This function writes a CBModel to file
    
     - *fba* an instance of an PySCeSCBM model
     - *fname* the output filename
     - *Dir* [default=None] use directory if not None
     - *separator* [default=','] the column separator
     - *only_exchange* [default=False] only output fluxes labelled as exchange reactions
     
    """    
    if Dir != None:
        assert os.path.exists(Dir), '\nPath does not exist'
        fname = os.path.join(Dir, fname)
    fname_r = fname + '.modelinfo.csv'
    objName = ''
    if len(fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions()) > 1:
        for J in fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions():
            objName += '%s_' % J
        objName = objName[:-1]
    else:
        objName = fba.objectives[fba.activeObjIdx].getFluxObjectiveReactions()[0]
    try:
        F = file(fname_r, 'w')
    except IOError:
        print '\nCSV file \"%s\" is locked by an external application (probably Excel) please close file and try again (or use a different filename).' % fname_r
        return 
    cntr = 0
    F.write('%s%s%s%s%s%s\n' % ('ObjectiveFunction',separator,objName,separator,separator,separator))
    F.write('\"%s\"%s%s%s%s%s\"%s\"%s\"%s\"%s\"%s\"\n' % ('Reaction',separator,'LowerBound',separator,'UpperBound',separator,'Name',separator,'Equation',separator,'Gene association'))
    for r in fba.reactions:
        if not only_exchange:
            GO = True
        elif only_exchange and r.is_exchange:
            GO = True
        if GO:
            if r.annotation.has_key('GENE ASSOCIATION'):
                gene = r.annotation['GENE ASSOCIATION']
            else:
                gene = 'none'
            bnds = fba.getReactionBounds(r.getPid())
            Lbnd = -numpy.inf
            Ubnd = numpy.inf
            if bnds != None:
                if bnds[1] != None:
                    Lbnd = bnds[1]
                if bnds[2] != None:
                    Ubnd = bnds[2]
            if r.reversible:
                equation = ' <==> '
            else:
                equation = ' >> '
            subs = ''
            prods = ''
            for rr in r.reagents:
                if rr.coefficient > 0.0:
                    if abs(rr.coefficient) == 1.0:
                        prods += ' + %s' % (rr.species_ref)
                    else:
                        prods += ' + %s %s' % (abs(rr.coefficient), rr.species_ref)
                else:
                    if abs(rr.coefficient) == 1.0:
                        subs += ' + %s' % (rr.species_ref)
                    else:
                        subs += ' + %s %s' % (abs(rr.coefficient), rr.species_ref)
            subs = subs[3:]
            prods = prods[3:]
            equation = subs + equation + prods
            F.write('\"%s\"%s%s%s%s%s\"%s\"%s\"%s\"%s\"%s\"\n' % (r.getPid(),separator,Lbnd,separator,Ubnd,separator,r.name,separator,equation,separator,gene))
            GO = False
    F.flush()
    F.close()
    print 'Reactions exported as CSV to %s' % fname_r


def printFBASolution(fba, include_all=False):
    """
    Prints the FBA optimal solution to the screen.
    
     - *fba* an FBA model object
     - *include_all* include all variables
     
    """
    OFflux = fba.objectives[fba.activeObjIdx].fluxObjectives[0].reaction
    OFvalue = fba.objectives[fba.activeObjIdx].value
    OFSense = fba.objectives[fba.activeObjIdx].operation
    print '\n\n**********\nModel: %s\n\n' % fba.getPid()
    print '%s objective: %s\nOptimal value: %s\n\n' % (OFSense,OFflux,OFvalue)
    if include_all:
        for J in fba.reactions:
            print '%s: %s' % (J.getPid(), J.value)
    print '**********\n'

def exportModel(fba, fname=None, fmt='lp', work_dir=None):
    """
    Export the FBA model in different formats:
    
     - *fba* the FBA model
     - *fname* [default=None] the exported filename if None then `fba.getPid()` is used
     - *fmt* [default='lp'] the export format can be one of: 'lp' (CPLEX), 'hformat' (Polyhedra), 'all' (both)

    Note that 'hformat' ignores 'fname' and only uses fba.getPid() this is a legacy behaviour
     
    """
    if work_dir == None:
        work_dir = os.getcwd()
    if fmt == 'all' or fmt == 'lp':
        WriteModelLP(fba, work_dir=work_dir, fname=fname)
    if fmt == 'all' or fmt == 'hformat':
        if fname == None:
            fname = fba.getPid().replace('.xml', '')
        WriteModelHFormatFBA2(fba, fname=fname, work_dir=work_dir)
        WriteModelHFormatFBA2(fba, fname=fname, work_dir=work_dir, use_rational=True)

def writeProteinCostToCSV(fba, fname):
    """
    Writes the protein costs 'CBM_PEPTIDE_COST' annotation toa csv file.
    
     - *fba* an instantiated FBA object
     - *fname* the exported file name
     
    """
    F = file(fname+'.csv','w')
    for R in fba.reactions:
        rid = R.getPid()
        pcost = ''
        avg_l = ''
        if R.annotation.has_key('CBM_PEPTIDE_COST'):
            pcost = R.annotation['CBM_PEPTIDE_COST']
        if R.annotation.has_key('CBM_PEPTIDE_LENGTH_MAX'):
            if R.annotation['CBM_PEPTIDE_LENGTH_MAX'] == None:
                avg_l = 1
            else:
                if R.annotation.has_key('CBM_AVG_PEPTIDE_LENGTH'):
                    avg_l = R.annotation['CBM_AVG_PEPTIDE_LENGTH']
                else:
                    avg_l = 0
        F.write('%s,%s,%s\n' % (rid,pcost,avg_l))
    F.flush()
    F.close()
    print 'Protein costs written to file: %s' % (fname+'.csv')

def WriteFVAtoCSV(id, fva, names, Dir=None, fbaObj=None):
    """
    Takes the resuls of a FluxVariabilityAnalysis method and writes it to a nice
    csv file. Note this method replaces the glpk/cplx_WriteFVAtoCSV methods.
    
     - *id* filename_base for the CSV output
     - *fva* FluxVariabilityAnalysis() OUTPUT_ARRAY
     - *names* FluxVariabilityAnalysis() OUTPUT_NAMES
     - *Dir* [default=None] if set the output directory for the csv files
     - *fbaObj* [default=None] if supplied adds extra model information into the output tables
     
    """
    if Dir != None:
        Dir = os.path.join(Dir, id+'.fva.csv')
    else:
        Dir = id+'.fva.csv'
    F = file(Dir, 'w')
    if fbaObj == None:
        F.write('name,optval,min,max,diff,red cost,minstat,maxstat\n')
    else:
        F.write('name,optval,min,max,diff,red cost,minstat,maxstat,"equation","subsystem","gene association","confidence level"\n')
    for Jidx in range(len(names)):
        if names[Jidx] != None:
            name = names[Jidx]
            optval = fva[Jidx][0]
            rc = fva[Jidx][1]
            min = fva[Jidx][2]
            max = fva[Jidx][3]
            diff = fva[Jidx][4]
            minstat = fva[Jidx][5]
            maxstat = fva[Jidx][6]
            if fbaObj == None:
                F.write('%s,%s,%s,%s,%s,%s,%s,%s\n' % (name, optval, min, max, diff, rc, minstat, maxstat))
            else:
                xInf = []
                Ro = fbaObj.getReaction(name)
                for k in ['Equation','SUBSYSTEM','GENE ASSOCIATION','Confidence Level']:
                    if Ro.annotation.has_key(k):
                        xInf.append(Ro.annotation[k])
                    else:
                        xInf.append('')
                F.write('%s,%s,%s,%s,%s,%s,%s,%s,"%s","%s","%s",%s\n' % (name, optval, min, max, diff, rc, minstat, maxstat, xInf[0], xInf[1], xInf[2], xInf[3]))

    F.flush()
    F.close()
    print 'FVA results written to: %s' % Dir


def WriteFVAdata(fva, names, fname, work_dir=None, roundec=None, scale_min=False, appendfile=False, info=None):
    """
    Takes the resuls of a FluxVariabilityAnalysis method and writes it to a nice
    csv file. Note this method replaces the glpk/cplx_WriteFVAtoCSV methods. Data is output as a csv file
    with columns: FluxName, FVA_MIN, FVA_MAX, OPT_VAL, SPAN
    
     - *fva* FluxVariabilityAnalysis() FVA OUTPUT_ARRAY
     - *names* FluxVariabilityAnalysis() FVA OUTPUT_NAMES
     - *fname* filename_base for the CSV output
     - *work_dir* [default=None] if set the output directory for the csv files
     - *roundec* [default=None] an integer indicating at which decimal to round off output. Default is no rounding.
     - *scale_min* [default=False] normalise each flux such that that FVA_MIN = 0.0 
     - *appendfile* [default=False] instead of opening a new file try and append the data
     - *info* [default=None] a string added to the results as an extra column, useful with `appendfile`
     
    """
    if work_dir != None:
        work_dir = os.path.join(work_dir, fname+'.fvadata.csv')
    else:
        work_dir = fname+'.fvadata.csv'
    if not appendfile:
        F = file(work_dir, 'w')
    else:
        F = file(work_dir, 'a')
    if info == None:
        F.write('%s,%s,%s,%s,%s\n' % ('Jid', 'min', 'max', 'optval', 'span'))
    else:
        F.write('%s,%s,%s,%s,%s,%s\n' % ('Jid', 'min', 'max', 'optval', 'span', 'info'))
    for Jidx in range(len(names)):
        if names[Jidx] != None:
            name = names[Jidx]
            max = fva[Jidx][3]
            min = fva[Jidx][2]
            optval = fva[Jidx][0]
            if roundec != None:
                max = round(max, roundec)
                min = round(min, roundec)
                optval = round(optval, roundec)
            if scale_min:
                if min > 0.0:
                    max = max-min
                    optval = optval-min
                    min = 0.0
                elif min <= 0.0:
                    max = max + abs(min)
                    optval = optval + abs(min)
                    min = 0.0
            if info == None:
                if roundec == None:
                    F.write('%s,%s,%s,%s,%s\n' % (name, min, max, optval, abs(fva[Jidx][3] - fva[Jidx][2])))
                else:
                    F.write('%s,%s,%s,%s,%s\n' % (name, min, max, optval, round(abs(fva[Jidx][3] - fva[Jidx][2]), roundec)))
                    
            else:
                if roundec == None:
                    F.write('%s,%s,%s,%s,%s,%s\n' % (name, min, max, optval, abs(fva[Jidx][3] - fva[Jidx][2]), info))
                else:
                    F.write('%s,%s,%s,%s,%s,%s\n' % (name, min, max, optval, round(abs(fva[Jidx][3] - fva[Jidx][2]), roundec), info))
    F.flush()
    F.close()
    print 'FVAdata results written to: %s' % work_dir

def writeSolutions(fname, sols=[], sep=',', extra_output=None):
    """
    Write 2 or more solutions where a solution is a dictionary of flux:value pairs:
     
     - *fname* the export filename
     - *sols* a list of dictionaries containing flux:value pairs (e.g. output by mod.getReactionValues())
     - *sep* [default=','] the column separator
     - *extra_output* [default=None] add detailed information to output e.g. reaction names by giving a CBModel object as an argument to *extra_output*.

    """
    assert len(sols) >= 2, "\nThere must be two or more solutions to work with"
    
    reac_ids = set([])
    for s in sols:
        reac_ids = reac_ids.union(set(s.keys()))
    reac_ids = list(reac_ids)
    reac_ids.sort()
    
    reac_names = []
    ##  reac_bnds = []
    ##  reac_eqns = []
    if extra_output != None:
        for r in reac_ids:
            reac_names.append(fba.getReaction(r).getName())
            ##  reac_bnds.append(fba.getReactionBounds(r.getPid()))

    F = file(fname+'.csv', 'w')
    for r in range(len(reac_ids)):
        row = '%s%s' % (reac_ids[r], sep)
        s_str = ''
        for s in range(len(sols)):
            if sols[s].has_key(reac_ids[r]):
                s_str += '%s%s' % (sols[s][reac_ids[r]], sep)
            else:
                s_str += '%s%s' % ('\"none\"', sep)
        row += s_str
        if extra_output:
            row += '\"%s\"\n' % reac_names[r]
        else:
            row = row[:-1]+'\n'
        F.write(row)
    F.flush()
    F.close()  
    print('\nSolutions written to: \"%s\"\n' % fname)    



